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Abstract. Computation of singular value decomposition (SVD) has been 
a topic of concern by many numerical linear algebra researchers. Fast 
SVD has been a very effective tool in computer vision in a number of 
aspects, such as: face recognition, eye tracking etc. At the present state 
of the art fast and fixed-point power efficient SVD algorithm needs to 
be developed for real-time embedded computing. The work in this pa¬ 
per is the genesis of an attempt to build an on-board real-time face and 
eye tracking system for human drivers to detect loss of attention due 
to drowsiness or fatigue. A major function of this on-board system is 
quick customization. This is carried out when a new driver comes in. 
The face and eye images are recorded while instructing the driver for 
making specific poses. The eigen faces and eigen eyes are generated at 
several resolution levels and stored in the on-board computer. The dis¬ 
criminating eigen space of face and eyes are determined and stored in 
the on-board flash memory for detection and tracking of face and eyes 
and classification of eyes (open or closed) as well. Therefore, fast SVD 
of image covariance matrix at various levels of resolution needs to be 
carried out to generate the eigen database. As a preliminary step, we 
review the existing symmetric SVD algorithms and evaluate their feasi¬ 
bility for such an application. In this article, we compare the performance 
of (1) Jacobi’s, (2) Hestenes’, (3) Golub-Kahan, (4) Tridiagonalization 
and Symmetric QR iteration and (5) Tridiagonalization and Divide and 
Conquer algorithms. A case study has been demonstrated as an example. 

Keywords: Fast SVD, Eigen Space, Face and Eye Detection 


1 Introduction 

SVD is a powerful tool in digital signal and image processing. It states that a 
matrix can be decomposed as follows, 

A = UEV^ (1) 

Where A^xn is a dense matrix, Umxm and Vnxn are orthogonal (or unitary) 
matrices and their columns are called left and right singular vectors respectively. 
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A" is a diagonal matrix and contains all singular values along its diagonal in a 
non-increasing order. For a symmetric matrix m = n and U and V span the 
same vector space. Hence computation of either U or V is sufficient. For any 
dense symmetric matrix Anxn, Eigen Value Decomposition (EVD) is defined as 
follows 


A = XAX'^ (2) 

Where X is the matrix of eigenvectors and A is the diagonal matrix containing 
eigenvalues along its diagonal. Eor symmetric matrices, eigenvalue decomposi¬ 
tions and singular value decompositions are closely related as follows 
Suppose that A is a symmetric matrix, with eigenvalues Ai and orthonormal 
eigen vectors Ui so that A = UAU'^ is an eigenvalue decomposition of A, with 
A = diag[XiX 2 ...Xn], U = [uiU 2 ---Un\ and UU'^ = I. Then an SVD of symmetric 
matrix A is, A = USV'^, where diagonal elements of S i.e. ai = abs{Xi) and 
Vi = sign{ai).Ui where sign(0) = 1. For symmetric positive definite matrices 
eigenvalue decomposition (EVD) and SVD leads to the same decomposition. 
Hence we will use eigenvalues/eigenvectors and singular values/singular vectors 
interchangeably. 

However, SVD has been an offline tool for digital signal and image process¬ 
ing applications for decades because of the computation complexity and memory 
requirement. Due to the increased resources of some of the recently introduced 
workstations, there are attempts to develop faster versions of SVD algorithm 
for real-time signal and image processing applications. The implementation of 
SVD in embedded platforms like DSPs, ARM and EPGAs is necessary for facil¬ 
itating efficient real-time image processing. Most of these platforms either have 
fixed-point processors or CLBs (Configurable Logic Blocks) to make the system 
cheaper and power efficient. Hence fast and fixed-point SVD algorithms are to 
be developed for such applications. The purpose of this work is to evaluate the 
existing SVD algorithms for their suitability on embedded platforms. 

In pattern recognition, eigenspace based method has been proposed for face 
tracking or face recognition m-m- To find the eigenspace, SVD (or eigenvalue 
decomposition) is used. There are several algorithms for SVD as stated in liter¬ 
ature M-m- Jacobi’s algorithm is known to be the oldest and slowest algorithm 
E-m .For symmetric matrices, though Jacobi’s algorithm generates accurate sin¬ 
gular values and singular vectors, the time of execution increases with the di¬ 
mension of the matrices and is only suitable as an offline tool.Two-sided and 
one-sided variants for Jacobi’s algorithm are stated in literature. Hestenes’ al¬ 
gorithm is a variant of one-sided Jacobi’s algorithm and is discussed in m-m- 
Being a one-sided version the time of computation is lesser than two-sided Ja¬ 
cobi’s algorithm. However, as the iteration is applied on the whole process this 
algorithm is also not suitable for online applications. Golub and Kahan proposed 
a two-step algorithm for computation of SVD. In the first step, a 

dense symmetric matrix is converted to a bidiagonal matrix, which is eventually 
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converted to a diagonal matrix using implicit QR iteration in the second step. 
Because the second phase of Golub-kahan algorithm is iterative in nature, it is 
much faster than Jacobi’s or Hestenes’ algorithm. A similar two step algorithm 
has been proposed for SVD, where a dense symmetric matrix is reduced to tridi¬ 
agonal matrix and then an implicit symmetric QR iteration is applied to reduce 
the symmetric tridiagonal matrix into a diagonal matrix. This algorithm is found 
to be faster and competitive with Golub-Kahan algorithm when a combination 
of QR and QL algorithm is used E-m- Still these algorithms could not satisfy 
real-time constraints as required by the signal and image processing platforms. 
Hence a Divide and Gonquer algorithm was proposed by J. J.M. Cuppen based on 
a rank-one modification by Bunch, Nielsen and Sorensen. This is the fastest algo¬ 
rithm till date when a complete eigensytem of a symmetric tridiagonal matrix is 
required [5] . A variant of the said algorithm by Gu and Eisenstat has been imple¬ 
mented in LAPAGK routine for matrices with dimension larger than 25 |13)-|14]. 

Faster performance is achieved when floating-point SVD is converted to fixed- 
point format and implemented in fixed-point platform [22) . Fast and Fixed-point 
SVD algorithm is also useful for reducing silicon area and power consumption 
in embedded platforms [25-[S5]. For Digital Signal Processing applications, at¬ 
tempts have been made to implement SVD algorithm using multiprocessor ar¬ 
rays [35] and GORDIC (GOordinate Rotation Digital Gomputer) based recon- 
figurable systems m-m- 


2 Existing Algorithms for SVD of Symmetric Matrices 

In this section, we analyse different SVD algorithms for symmetric matrices 
along with their complexity, advantages and disadvantages. 

1. Jacobi’s Algorithm |5] 

2. Hestenes’ Algorithm [10] 

3. Golub-Kahan Algorithm [TT|-[T3] 

4. Tridiagonalization and Symmetric QR iteration E-m 

5. Tridiagonalization and Divide and Gonquer Algorithm [^ 

6. Bisection and Inverse Iteration [S] 


2.1 Jacobi’s Algorithm 

Jacobi’s algorithm is the oldest and slowest available method for computing SVD 
by implicitly applying iteration on a dense symmetric matrix A. The method is 
more or less similar to the method for eigenvalue decomposition of symmetric 
matrices. Jacobi’s method computes SVD of a symmetric matrix A with higher 
relative accuracy when A can be written in the form A=DX, where D is diagonal 
and X is well conditioned [5]. Thus 

J^AJ = S 


(3) 
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Table 1. Computational complexities of SVD algorithms 


Algorithms 

Time Complexity 

Jacobi’s Algorithm 


Hestenes Method 

0{n^) 

Golub-Kahan Algorithm 

(V[(|n'’’ -1- 0{n‘‘))] and 

(Bidiagonalization + Implicit QR Iteration) 

U,V[An^ + 0{'tA)]) 

Tridiagonalization + Symmetric QR Iteration 

9,W + 0{n^) 

Tridiagonalization + Divide and Conqner Method 

|n® -I- 0(n^) 


In each step we compute a Jacobi rotation with J and update A to AJ, where 
J is chosen in such a way that two off-diagonal entries of a 2 x 2 matrix of A 
are set to zero in j"’'AJ . This is called two-sided or classical Jacobi method. 
Again, this can be achieved by forming G = A"^ A and performing iteration over 
G instead of A. The eigenvalues of G are the square of the singular values of A. 
Since J^GJ = J^A'^AJ = (AJ)^(AJ), we can obtain E by merely computing 
AJ . This is termed as one-sided Jacobi rotation. Though Jacobi’s algorithm for 
calculating SVD is the slowest among all presently available algorithms, it pro¬ 
duces relatively accurate singular values and singular vectors. Jacobi’s algorithm 
can be implemented in parallel as the individual steps are not interdependent. 
Parallel Jacobi’s algorithm is discussed in [15]. 


Table 2. Two-sided Jacobi’s algorithm for SVD 


Q — I (/ = identity matrix) 
repeat 

s=c.t 

for J = 1 : (n — 1) 

= c,J{k,k) = c, J{j,k) = s,J{k,j) = -s 

for k = {j + 1) ■. n 

A = J'^AJ 

if \A{j, A:)| is not too small 

Q = Q.J 

r _ (A(k,k)-A(jJ)) 

S 2A(j,k) 

endif 

if 5 = 0 

endfor 

^ (ien-v/i+?2' 

endfor 

else 

ak = \A{k, k)\ 

^ _ signd) 

U = Q 

endif 
c = . ^ 

Vk = sign{A{k, k)).Uk 

V = [viV 2 ■ ■ - Vn] 

y/l + t2 

contd ... 



2.2 Hestenes’ Algorithm 

Hestenes’ algorithm is based on one-sided Jacobi’s algorithm. A symmetric ma¬ 
trix of size n X n is used to generate an orthogonal rotation matrix V so 
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that the transformed matrix A' = AV = W has orthogonal columns. Now if 
we normalize each non-null column of matrix W to unity, we get the relation 
W = AV = US = U■diag{ai,a 2 , cr„), A = USV'^, V = 0"=! singular 

values of A are ai = || a' II 2 , where A' = ... a'^. 


Table 3. Hestenes’ algorithm for SVD 


V = I 


repeat until convergence 

c = —J=, s = t.c 

y/l+t? 

for p = 1 : (n — 1 ) 

end 

lor q = {p + 1) ■. n 

J=I 

a = Ai:,pfAi:,p),d = Ai:,qfA{-.,q) 

J{p,p) = c, J{q,q) = c,J(p,q) = s,J{q,p) = -s 

7 = A{:,p)'^ A{:,q) 

A = A.J, V = V.J 

if 7 = 0 

endfor 

c = 1 ; s = 0 

endfor 

else 

for A: = 1 ; n 

II 

ai = P(:,fe)|| 

^ sign(C) 

end 

dCI+Vi+d") 


contd ... 

U = AVE-^ 


2.3 Golub-Kahan Algorithm 

This algorithm can be segmented in two phases, 

i) In the first phase a dense symmetric matrix is converted to a bi-diagonal ma¬ 
trix by orthogonal similarity transformation. 

ii) This bi-diagonal matrix is then converted to diagonal matrix using Implicit 
QR iteration. 

The standard algorithm may compute singular values with poor relative ac¬ 
curacy. However, with modified algorithm by Demmel-Kahan smaller singular 
values may be computed with high relative accuracy m- This algorithm has 
two steps as shown in equation 


/* =1= * 


/* * 0 o\ 


/* 0 0 0\ 

* * * * 

UiAVU=B^ 

0**0 

U2BV2'^ = S^ 

0*00 

* * * * 

STEP I 

0 0 * * 

STEP II 

0 0*0 



P 0 0 *^ 


poo*/ 


Hence, U 2 ^BV 2 = U 2 ^Ui'^AViV 2 = {UiU 2 f AV 1 V 2 = U'^AV = E and S is 
the diagonal matrix containing singular values of the symmetric matrix A. In bi- 
diagonalisation process we have used householder method to make the elements 
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Table 4. Golub-Kahan algorithm for SVD 


function [u, (T]=houszero(x) 
m = max{\xi\), i = 1,2, ... ,n 

Ui = — ,i = 1,2, ... ,n 
sign{0) = 1 

cr = sign{ui)\/ui^ + U 2 ^ + • • • + Un^ 
ui = ui + a 

a — —m.a 

end houszero 

function [c, s, r]=rot(f,g) 

if / = 0 then c = O', s = V, r — g 

elseif (I/I > \g\) then 

t = / i P = -I- c = = t.c; r = tt.f 

else 

t = ^-,11 = Vl + H',s= -^',c = t.s; r = tt.g 

endif 

end rot 

Golub-Kahan algorithm for SVD 


Ui=I,Vi=I 

repeat 

for fc = 1 : (n — 1) 

for i = 1 : (n — 1) 

[u, cr] = houszero(A(fc : m, k) 

U2=I 

HI = 1-2^ 

[c, s, r] =rot{A{i, i),A{i, i + 1)) 

Pi=I 

Q = I 

P\{k ■. m,k ■. n) — HI 

Q{i : i + l,i : i + 1) = ( '^ * ) 


V-s cj 

A{k : m,k : n) — Hl.A{k : m,k : n) 

A = AQ'^ 

Ui=Ui.Pi 

V2 = V2.Q'^ 

if fc < (n — 2) 

[c, s, r] =rot{A{i, i),A{i + 1, i)) 

[u, a] =houszero(A(fc, fc + 1 : n)^) 

Q^I 

/ \ 

H2 = 7 - 

v-^ V 

Q{i:i + l,i:i + l) ^ 

P2 =7 

A = QA 

P 2 (fc : m — 1, fc : n — 1) = 772 

U2 = U2.Q 

A{k : m,k + 1 : n) — A{k : m, fc + 1 : n).H2 

endfor 

Vl=Vl.P2 

E = abs{A) 

endif 

= U 1 .U 2 

endfor 

V = V 1 .V 2 

contd ... 
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of a column or row zero. Another optimized variant of this algorithm is also 
available and is known as Golub-Kahan-Chan algorithm. 

2.4 Tridiagonalization and Symmetric QR iteration 

In this method, the dense symmetric matrix is first converted to symmetric tridi¬ 
agonal matrix in finite number of steps and the tridiagonal matrix is eventually 
converted to diagonal form by symmetric QR iteration. When a combination of 
QL and QR algorithm instead of QR is used, a stable variant is obtained. The 
steps of this algorithm is stated in equation ([5|) . 

Method of tridiagonalization 

Householder reflection can be used for tridiagonalization. In this method, one 
row or column is picked up, and householder matrix is reconstructed to make all 
the elements of the row or column zero except the first one [5]. 

Hessenberg reduction m] is a process by which a dense matrix is converted 
to an upper or lower Hessenberg matrix by orthogonal similarity transforma¬ 
tion and thus does not change the eigenvalues or singular values. For a dense 
symmetric matrix, this Hessenberg reduction produces an upper as well as a 
lower Hessenberg matrix or a symmetric tridiagonal matrix which is used in 
symmetric QR iteration process to produce a diagonal matrix.In this process, 
the eigenvalues or singular values remain same as that of the original dense sym¬ 
metric matrix. Householder or Givens method may be applied for Hessenberg 
reduction. We have used Givens method to produce symmetric tridiagonal ma¬ 
trix. 

Lanczos method is useful to transform a symmetric dense matrix to a sym¬ 
metric tridiagonal matrix. This method suffers from loss of orthogonality among 
Lanczos vectors with increased number of iterations. For this reason reorthogo- 
nalization is required to obtain proper orthogonal vectors. 

Symmetric QL and QR iteration : 

This is an iterative process and has a complexity of 0{n^). After reducing the 
dense symmetric matrix into its tridiagonal form by orthogonal similarity trans¬ 
formation, symmetric QL or QR iteration is performed depending on which of 
the first and last diagonal elements of the symmetric tridiagonal matrix is larger. 
If the first diagonal entry is larger than the last one, QR is called upon to per¬ 
form, else QL is called. 


A * =1= *\ 


/* * 0 o\ 


A 0 0 o\ 

* * * * 

UiAVi'^=T^ 

* * ^ 0 

U2TV2'^ = S^ 

0*00 

* * * * 

STEP I 

0 * * * 

STEP II 

00*0 

\****) 


\0 0 * 


\poo*J 
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Table 5. Tridiagonalization and Symmetric QR Iteration 


function [c, s] =givens(a, 6) 

Pi{k + 1 : n, k + 1 : n) = H 

if 6 = 0 

A = PiAPi 

c = 1; s = 0 

Pi = Pi.Pi 

else 

end 

if |6| > |a| 

repeat 

r=-T 

j C^n-l.n-l—in,n) 

2 

1 



^ (d+sig-a(d)^d^+tl^^_^) 

else 

X = til — fj.-, Z = t21 

^ = -f 

for fc = 1 : (n — 1) 

c = , ^ ; s = - 

[c, s]=givens(a;, z) 

endif 

T = {Gl)TGk, where Gk = G{k, k + l,e) 


Pi = PiGfc 

Ui=I 

if fc < (n — 1) 

for fc = 1 ; (n — 2) 

X = tk-\-l,k 

[u, a] =houszero(a;) 

z = tk+ 2 ,k 

H - h 

{u-'- u) 

endif 

Pi =7 

endfor 

contd.. 

P = abs{T),U = Pi 


2.5 Tridiagonalization and Divide and Conquer Algorithm 

The method was first proposed by J.J.M. Cuppen [14] based on a rank-one mod¬ 
ification by Bunch, Nielsen and Sorensen [18] . However, the algorithm became 
popular after a stable variant for finding singular vectors or eigenvectors was 
found out in 1990 by Gu and Eisenstat This algorithm is the fastest SVD 
algorithm available till date. A dense symmetric matrix is first converted to sym¬ 
metric tridiagonal matrix Then the symmetric tridiagonal matrix is 

divided into two parts by rank one update and again each of the smaller matrices 
is divided till sufficiently smaller (matrix dimension = 25) matrices are formed. 
Then QR and QL iteration may be applied to find the SVD of the smaller matri¬ 
ces and using rank one update smaller solutions are combined together to form 
the complete SVD of the symmetric tridiagonal matrix. With a combination of 
previously stated steps, a complete eigensystem of a dense symmetric matrix 
can be found out. 

Two major parts for finding the eigensystem of a symmetric tridiagonal ma¬ 
trix are divide and conquer. It works by breaking down a problem into two or 
more subproblems of the same type until the subproblem becomes simple enough 
to be solved directly. The solutions to the subproblems are then combined to gen¬ 
erate a solution to the original problem. The most significant part is the solution 
of secular equation. The solution of the secular equation involves function ap¬ 
proximation for finding the desired roots. For matrices with dimension greater 
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than 25 this is the fastest method for finding the complete eigensystem of a 
symmetric tridiagonal matrix till date [S] . 


Structure of the Algorithm [5] This includes dividing the symmetric tridi¬ 
agonal matrix into two parts by removing the subdiagonal entry by rank one 
modification. Now with the known eigensystem of the two new symmetric tridi¬ 
agonal matrices the secular equation is constructed. Solution of secular equation 
gives the eigenvalues of the original matrix from which the eigen vectors can also 
be computed. A dense symmetric matrix is converted to a symmetric tridiagonal 
matrix and then the divide and conquer algorithm is applied as m- 

1. Step 1 - Dividing the symmetric tridiagonal matrix into smaller matrices by 
rank one modification. 

2. Step 2 - Sorting eigenvalues and eigenvectors obtained from smaller matrices 
in an increasing order using permutation matrix. 

3. Step 3 - Deflation (when updation of eigenvalues and eigen vectors is not 
required) due to smaller coefficient or smaller components in the vector used 
for rank one modification and also for repeated eigenvalues are considered. 

4. Step 4 - Formation and solution of secular equation using non-deflated eigen¬ 
values. 

5. Step 5 - Combining the solution obtained from secular equation solver and 
deflated eigenvalues to obtain the complete eigensystem of the symmetric 
tridiagonal matrix. 

Forming the symmetric tridiagonal matrix from a dense symmetric matrix 

U'^AU = T (6) 

A - Dense symmetric matrix 
T - Symmetric tridiagonal matrix 

U - Eigenvector matrix obtained from tridiagonalization process 


Divide and Conquer Algorithm 

Step 1 Rank one modification of the symmetric tridiagonal matrix 


T = 


Ti 0 

0 T2 


T 

■iVV = 


QiDiQi^ 

0 


T = 


Qi 0 

0 Q2 


Di 0 

0 D2 


0 

< 52 ^ 2 ( 52 ^ 


Qi 

0 


iVV 


Q2 


(7) 

( 8 ) 


Where, T\,T 2 are smaller tridiagonal matrices. 

D is a diagonal matrix containing the eigenvalues of the smaller tridiagonal 
matrices. 

Qi,Q 2 are eigenvector matrices after eigendecomposition of smaller tridiagonal 
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matrices. 

_ f Qi^ 0 \ / last column of Qi^\ 

^ V 0 ^ \first column of Q 2 ^) 

Step 2 Sorting using permutation matrix 

Following equation is used for sorting process to arrange the eigenvalues and 
eigenvectors in an increasing order. 

D + = P^P{D + puu^)P^P (9) 

Where P is the permutation matrix. 

Step 3 Reducing the computation using deflation 
Deflation occurs when 

1 . p is very small, 2 . smaller weights (ui) due to smaller components in the vector 
used for rank one modification and 3. multiple eigenvalues m- When deflation 
occurs eigenvalues and eigenvectors do not require updation. Hence a saving in 
computation is achieved. 

Step 4 Formation and solution of secular equation 
\{{D + puu^) - XI)\ = 0 
\{{D - XI){I + p{D - XI)-^uu'^)\ = 0 
Since, \{D — A/)| 7 ^ 0, |(/ + p{D — XI)~^uu^)\ = 0 
Now, 

n 2 

|((/ + p{D - XI)~^uu'^)\ = 1 + pu^{D - XI)~^u = 1 + 

t—1 

^ + bmuu"’' = [i?] [A] [i?^] [obtainedby solving secular equation) 

Eigenvalues and eigenvectors are obtained by solving a secular equation like 
equation (nni) [I3-II1]- 

Step 5 Combining the solutions from Step IV and V the complete eigensystem 
of the symmetric tridiagonal matrix is obtained 

^ {{^0 D,) + = {QmiMRflQf = ( 12 ) 

Where, [X] = [Q][i?]. 

Hence, 

A = [U][T][uf = [U][X][A][xf[uf = [H][H][H]^ (13) 

There are some issues related to divide and conquer algorithm. They are dis¬ 
cussed below. 

Sorting with Permutation Matrices j20j If di < d 2 ■ ■ ■ < dn then the 
sequence of eigenvalues obtained will be Ai < A 2 < • • • < A„. However, we may 
not come across a diagonal matrix with sorted diagonal elements after rank one 
modification and QR or QL iteration. Hence, we need to apply permutation to 
sort them in ascending order using permutation matrix. 


fDi 0 
\0 D 2 


( 10 ) 

( 11 ) 
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Let, D = 


Let, u = 


/13.1247 0 

0 201.9311 


0 0 

0 0 

0.0693 0 

0 26.7189/ 


0 

V 0 

^- 0 . 5421 ^^ 

-0.4540 

0.2128 

\-0.6743/ 


Now applying permutation using permutation matrix, P = 


/0.0693 0 


equation (ED , we obtain D (sorted) = 


13.1247 

0 

0 


u(modified)= 


f 0.2128 

-0.5421 

-0.6743 

\-0.4540/ 


/O 0 1 0\ 

10 0 0 
0 0 0 1 
\0 1 0 0 / 

0 0 

0 0 

16.7189 0 

0 201.9311/ 


and using 


and 


Solving Secular Equation The roots of the secular equation 1 + pY^^=i d -\ 
are the eigenvalues of the original matrix, p (rho) is the sub diagonal entry 
which creates the rank one modification. Ui^ is the weight over the pole. The 
secular equation has poles at the eigenvalues of D and zeros at the eigenvalues 
of H + puvT'. 

According to interlacing property : 

1. If rho is greater than zero, the roots lie in such a manner that : di < Ai < 
d2 < X 2 < ■ ■ ■ < dn < Xn 

2. If rho is less than zero then : Xi < di < X 2 < d 2 < ■ ■ ■ < Xn < dn 


Assuming rho is greater than zero for i < n, the roots lie in between di and 
di+i, but for i = n, the root lies in a manner that dn < Xn < dn + puu^ [T5] . 
"16.7118 10.7270\ 


For the given matrix 


^10.7270 34.2341/ 
zero, the nature of the roots are examined below 

The eigenvalues of the previously stated matrix are 


where rho (10.7270) is greater than 


11.6228 0 
0 39.3231 


where the diagonal matrix D is 


5.9848 0 

0 23.5071 


. In Figure [T]we can see bold 


lines and vertical dotted lines. The points, where the bold lines cross the real 
axis at zero, are the roots for the corresponding matrix. The vertical dotted lines 
represent the diagonal elements after the rank one modification and eigendecom- 
position of smaller matrices. 
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Fig. 1. Roots showing interlacing property with rho > 0 


For a matrix ^ ^io'^7270 34^2341^) (-10.7270) less than zero, the roots 

appear as shown in Figure [D 

The eigenvalues of the above said matrix with rho < 0 are 


11.6228 0 
0 39.3231 


where the diagonal matrix D is 


27.4388 0 

0 44.9611 


is clearly found from Figures [T]and [2] 


. The interlacing property 


Methods of Solving Secular Equation |18) We assume that p (rho) is 
greater than zero. The solution of secular equation is based on function approx¬ 
imation. So we may think of using Newton’s iterative procedure to solve the 
equation, but Newton’s method is based on local linear interpolation. At some 
points of initial guess (Aq) for the desired root, the linear approximation would be 
horizontal and the next approximation would be a large negative number which 
is not useful [5]. This happens particularly when the weights are small. Therefore, 
we go for rational osculatory interpolation. Rational osculatory interpolation of 
secular function is the combination of two kinds of rational functions. If the 
secular function is /(A), the two rational functions are respectively ‘^(A) and 
<F(A). ^(A) is the sum of positive terms and *F(A) is the sum of negative terms. 
The inequality among them is —oo < ’F(A) < 0 < ‘^(A) < -boo. In the following 
section, we will explore a set of schemes for solving secular equation. The first 
scheme is named as approaching from left because the algorithm will produce 
a sequence of monotonically increasing approximations to desired root provided 
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liiteilncing pio|ieity|witli p< 0) 



Rootsofseculaifuiictioii(x| 


Fig. 2. Roots showing interlacing property with rho < 0 


the initial guess is less than the desired root. Starting with the initial guess which 
is less than the desired root, the scheme will yield a sequence of approximations 
converging monotonically upwards to the desired root. On the other hand, if 
the initial guess exceeds the required root so much that the next approximation 
does not converge, the scheme fails to get the desired root. Approaching from 
left scheme proves better for i < n, but for i = n according to the interlacing 
property the root lies in the range < A„ < Here the initial guess 

exceeds the desired root and the second scheme i.e. approaching from right will 
yield a sequence of approximations converging monotonically downwards to de¬ 
sired root. A middle way method has been adopted by combining both the above 
schemes. This scheme takes both nearby poles into consideration, while the other 
two schemes take each of the nearby poles into consideration. Here the scheme 
handles two parts. The first part is, if the combination of two rational functions 
is greater than zero, the root is closer to di and if the function is less than zero, 
the root is closer to Interpolation of both >f'(A) and d>{\) is carried out by 
taking both nearby poles into consideration. Following [5], 


h{X) = 



C2 

di+1 — A 


+ C3 


(14) 


/(A) = 1 + ELi + EL .+1 

= 1 -f ^(A) -f !F(A) (15) 

Where coefficients ci,C 2 and C 3 are computed by interpolating /(A). In [TB] 
the above mentioned secular function has been represented similarly using the 
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system’s matrix as: 


D H —uu 
P 


T 


/l(A) — p T + + 


+ 


S' 


dk — X dk+i — A 


(16) 

(17) 


Where r + is approximated to 'f'(A) and R + is approximated to 

<?(A). The initial guess (Aq) plays a very significant role in finding solution for 
secular equation. For k < n the initial guess is calculated as and for 

k = n the initial guess is d„ + puu^. So with Aq, we calculate the value of /i(A), 
and if the rational function is greater than zero, the root is closer to di and if it 
is less than zero, the root is closer to di+i. So according to the placement of the 
root, if it is close to di then we shift each of Aq, di and di+i by subtracting di 
from them. With the new value of the initial guess (Ai), we solve the quadratic 
equation formed by the two rational functions, and rj (iterative correction ) is 
found. The desired root is found by computing (Ai + rj). We know that the 
weights over the poles play a major role in the solution of secular equation. One 
of the circumstances is associated with weights. In order to overcome this, fixed 
weight method was implemented with the structure same as middle way method. 
By combining middle way method and fixed weight method, a hybrid scheme was 
designed to make iteration faster. In this scheme, the function interpolation was 
carried out with dk, dk+i and dk-i- The rational function for hybrid scheme is 
given below. 


h{X) = c + 


Uk 


dk-i — X dk — X dk+i — X 


(18) 


Where c = p + r + R. Once the eigenvalues are obtained, eigenvectors can be 
computed using the following equation 


{XI -D)-^u 
II {XI-D)-\ II 


(19) 


Unfortunately, the above equation does not produce accurate eigenvectors [17,20] 
and the following equation is used for modifying u vectors where u vectors are 
updated and used in equation m 


Uk 


Y\j=i{dk Aj) f([^._^.(Aj dk) 

\| PU^Zlidk - dj) U]=k+iid3 - dk) 


( 20 ) 


In this way, the eigenvectors computed (using equation [T51) are relatively accu¬ 
rate [T7] . 


Apart from the above SVD algorithms. Bisection with Inverse Iteration can 
also be used. This is specially used when singular values or eigenvalues are re¬ 
quired within a specified range. This algorithm suffers from loss of orthogo¬ 
nality among the computed singular vectors or eigenvectors. Hence forced re- 
orthogonalization is necessary. 














Fast SVD 


15 


3 Case Study : Real-time Face and Eye Tracking 

A fast SVD based face tracking scheme based on eigen space is shown in Figure |31 
The scheme is composed of two sections - hardware and software. Fast and fixed- 
point SVD is required in this scheme for rapid customization of the system for a 
particular human subject. There are two paths: path 1 is for training and path 
2 is for tracking purposes. The subject is asked to put its face within a box area 
(frame) in the image. The subject is then required to follow some instructions 
which direct him to make different poses with different face orientations. The 
training set is prepared online from extracted face images with different face 
orientations (front, front up, front down, looking left, looking right etc.) and it is 
customized for the person sitting in front of the camera. The covariance matrix is 
constructed from extracted face images with different poses at a frame rate of 30 
fps. Fast SVD is performed on the covariance matrix to prepare the eigen space or 
feature space of face or eye in real-time. Now face and eye detection and tracking 
is carried out by projecting a block image from the incoming frame (from web 
camera) onto the feature vector space and comparing the reconstructed image 
with a standard face or eye image. Face and eyes are detected at approximately 
10 fps in the laboratory environment. 



Fig. 3. Schematic diagram of the proposed online training based face and eye tracking 
system. 


The steps for the scheme mentioned above are given below: 

Eigenfaces are calculated as in [3]: 

Step - 1: Obtain face images (for training), each of dimension say 
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Fig. 4. EigenFaces (top), Eigen Eyes (for open and closed eyes) and Eye Detection 
using Fast SVD based algorithm (bottom). 
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N X N. 

Step - 2: Represent every image li as a vector Fi (of dimension x 1). 

Step - 3: Compute the average or mean face vector W: 

M 

‘'“mEC (21) 

i=l 

Step - 4: Subtract the mean face from each face vector Ff. 

^^=F,-F (22) 

Step - 5: The covariance matrix C given by: 

M 

C == AA^ {N^ X N^) (23) 

n—1 

where A= [<j)i^ (j> 2 , 4>m] x M matrix) is very large, compute A instead 

(M X M matrix, M ^ N). 

Now AA^ = USE^U’^ and A'^A = VS'^EV^. 

Step - 6: Compute the singular vectors Vi of A^ A. Using the equation Ui = Avi 
[14], singular vectors Ui of AA"'" are obtained. 

Step - 7: Depending on energy content of the components, keep only K singular 
vectors or eigenvectors corresponding to the K largest singular values. These K 
singular vectors are the eigenfaces corresponding to the set of M face images. 
Step - 8: Normalize these K eigenvectors so that Ui = 

Step - 9: Civen a test image F , subtract the mean image F : (p = F — F. 
Project the normalized test image onto the face space and obtain weight vector, 
C = U'^F, where 17 = [a;iW 2 • ■ • Wfcj- 

Step - 10: Compute reconstructed image S = WiUi, {wi = uA<l>) . 

Step-11: Compute error, e = WF — S\\. 

Step-12: Classify test image as belonging to face class image for which the error 
norm e is minimum. 

PERcentage eye CLOSure (PERCLOS) is defined as the proportion of time for 
which the eyelid remains closed more than 70-80% within a predefined time 
period. Level of drowsiness can be judged based on the PERCLOS threshold 
value [53]. Let Na be the number of eye frames belong to the open or attentive 
category out of Nm number of eye frames captured in a minute. Hence {Nm — 
Na) is the number of eye frames belonging to the inattentive category. Then 
PERCLOS value per minute is 

PERCLOS = ~ X 100% (24) 

However, the PERCLOS value computed per minute is not the correct measure¬ 
ment of drowsiness. Literature suggests that 20 minutes bout to bout measure¬ 
ment is a prominent indicator of drowsiness |29j , |30j . So a trade-off between ex¬ 
ecution time and accuracy of PERCLOS is requird. PERCLOS measured within 
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a 3 minute time interval is found to indicate the level of drowsiness reasonably 
well. A threshold over the PERCLOS value can be used to indicate level of 
drowsiness [30] . A comparison of eigenspace preparation using QR and QL and 
divide and conquer algorithm is shown in Table |3| 

Table 6. Comparison of Execution Time (sec) of QL and QR and Divide and Conquer 
Algorithm 


Order 

QR and QL Algorithm 

Divide and Conquer Algorithm 

50 

0.046875 

0 

100 

0.046875 

0.03125 

200 

0.15625 

0.03125 

500 

1.671875 

0.046875 

1000 

9.421875 

0.09375 


4 Conclusion and Future Scope 

In this article different symmetric SVD algorithms have been evaluated for real¬ 
time preparation of eigen space for face and eye tracking to assess the level 
of alertness of human driver. Different symmetric SVD algorithms have been 
discussed with their complexity, advantages and disadvantages. Issues related to 
fast Divide and Conquer SVD algorithm have been discussed. A case study of 
real-time face and eye tracking has been illustrated with comparison of execution 
time for eigen space preparation with QR and QL and Divide and Conquer SVD 
algorithms. The scheme has been implemented in a workstation with Intel quad 
processor (2.5 GHz) and 3 GB DDR RAM. 

The present work has been undertaken with an aim to implement the algorithm 
in an embedded platform where processing speed and other constraints exist. 
Implementation of fixed-point fast SVD for face and eye tracking could be useful 
for fixed-point DSPs and other fixed-point embedded platforms like ARM. To 
convert floating-point algorithms into fixed-point,it is necessary to first estimate 
the range of the required floating-point variables followed by choice of Q-format. 
Then the floating-point variables and floating-point arithmetic equations are 
converted to fixed-point format [11],[23]. The fixed-point algorithm is then tested 
for error and the algorithm is optimized for our purpose. This fast and fixed-point 
SVD algorithm could be used for other embedded signal and image processing 
applications.Development of other fixed-point signal processing algorithms like 
Wavelet Transform, Hidden Markov Model could be few future extension. 
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